setwd("~/Dropbox/research/accommodation/replicationfiles")

pacman::p_load(tidyverse,estimatr,margins,texreg,huxtable,ggeffects,interflex)

# Load data ----

acdf <- read_rds("accommodation-repdata.rds")

# Evaluate immigration index reliability ----

relstats <- acdf %>% 
  dplyr::select(att_imm1,att_imm2) %>% 
  as.data.frame() %>% 
  psych::alpha() 

relstats$total$std.alpha

# Load helper functions ----

source("helperfunctions.r")

# Manipulation check -----

## Fig 1 ----

# manipulation checks
mcm_biv <- lm(sdimmpos~expcondition,data = acdf)
mcm_adj <- lm(sdimmpos~expcondition+education+occupation+partygroup15,data = acdf)

plotexpmargins(mcm_biv,mcm_adj,"Effect on perceived SD restrictiveness","Manipulation check")

ggsave("fig1.pdf",width=5,height=3)


# Main results -----

## Fig 2a -----

# SD ptv
sdptvm_biv <- lm(sdptv~expcondition,data=acdf)
sdptvm_adj <- lm(sdptv~expcondition+education+occupation+partygroup15,data=acdf)

plotexpmargins(sdptvm_biv,sdptvm_adj,"Effect on Social Democrats PTV","Social Democrats propensity to vote")

ggsave("fig2a.pdf",width=4,height=4)


## Fig 2b ----

#for context: how often does max left bloc ptv differ from sd ptv?
ptvvars <- select(acdf,sdptv:nbptv) %>% 
  mutate(maxlbptv=pmax(sdptv,rvptv,sfptv,elptv,na.rm=T),
         maxrbptv=pmax(veptv,dfptv,nbptv,na.rm = T),
         lbpref=ifelse(maxlbptv>maxrbptv,1,0)) %>% 
  filter(lbpref==1) %>% 
  mutate(sdtopptv=ifelse(sdptv>=pmax(rvptv,sfptv,elptv,na.rm=T),1,0))

table(ptvvars$sdtopptv) %>% prop.table()
#i.e., among respondents whose highest PTV is for a left bloc party, 56 pct. have another left wing party as top PTV

#max left bloc ptv
maxlbptvm_biv <- lm(maxlbptv~expcondition,data=acdf)
maxlbptvm_adj <- lm(maxlbptv~expcondition+education+partygroup15,data=acdf)

huxreg(maxlbptvm_biv,maxlbptvm_adj)

plotexpmargins(maxlbptvm_biv,maxlbptvm_adj,"Effect on highest left bloc PTV","Highest left bloc propensity to vote")

ggsave("fig2b.pdf",width=4,height=4)

#comparing accommodation effect to adversarial
lm(maxlbptv~expcondition_refadv,data=acdf) %>% tidy()
lm(maxlbptv~expcondition_refadv+education+occupation+partygroup15,data=acdf) %>% tidy()

huxtable::huxreg(mcm_biv,sdptvm_biv,maxlbptvm_biv)

# Heterogeneous effects ----

## Fig 3a ----

aiintm_biv <- lm(sdptv~antiimm*expcondition,data=acdf)
aiintm_adj <- lm(sdptv~antiimm*expcondition+gender+agegrp+education+partygroup15,data=acdf)

pgintm_biv <- lm(sdptv~expcondition*partygroup15,data=acdf)
pgintm_adj <- lm(sdptv~expcondition*partygroup15+gender+agegrp+education,data=acdf)

aiintmdf <- bind_rows(getintmarginsdf(aiintm_biv,moderator="antiimm",refadv = F,modeltype="Bivariate",comptype="Accomodative v. Control"),
                    getintmarginsdf(aiintm_adj,moderator="antiimm",refadv = F,modeltype="Adjusted",comptype="Accomodative v. Control")) %>% 
  mutate(model=fct_relevel(model,"Bivariate"),
         comp=fct_rev(comp))

pgintmdf <- bind_rows(getintmarginsdf(pgintm_biv,moderator = "partygroup15",refadv = F,modeltype="Bivariate",comptype="Accomodative v. Control"),
                      getintmarginsdf(pgintm_adj,moderator = "partygroup15",refadv = F,modeltype="Adjusted",comptype="Accomodative v. Control")) %>% 
  mutate(model=fct_relevel(model,"Bivariate"),
         pg=fct_relevel(partygroup15,c("Left bloc","Other right")),
         comp=fct_rev(comp)) %>%
  filter(pg %in% c("Left bloc","Other right","DPP"))

#get plot yaxis limits
min(aiintmdf$lower)
max(pgintmdf$upper)

#plots
ggplot(data=aiintmdf,aes(x=antiimm,y=AME,group=model,color=model)) +
  geom_hline(yintercept = 0,linetype="dotted") +
  geom_point(position = position_dodge(width=.05)) +
  geom_rug(data=acdf,aes(x=antiimm,y=0),sides="b",alpha=.04,position="jitter",inherit.aes = F) +
  geom_errorbar(aes(ymin=AME-1.65*SE,ymax=AME+1.65*SE),width=0,size=1.2,position = position_dodge(width=.05)) +
  geom_errorbar(aes(ymin=AME-1.96*SE,ymax=AME+1.96*SE),width=0,size=.8,position = position_dodge(width=.05)) +
  theme_bw() +
  labs(x="Anti-immigration attitude",y="Effect on PTV") +
  ggtitle("Social Democrats propensity to vote") +
  ylim(c(-.09,.25)) +
  scale_color_manual(values=c("black","gray")) +
  theme(legend.position = "none")

ggsave("fig3al.pdf",width=6,height=4)

ggplot(pgintmdf,aes(x=pg,y=AME,group=model,color=model)) +
  geom_hline(yintercept = 0,linetype="dotted") +
  geom_point(position = position_dodge(width=.25)) +
  geom_errorbar(aes(ymin=AME-1.65*SE,ymax=AME+1.65*SE),width=0,size=1.2,position = position_dodge(width=.25)) +
  geom_errorbar(aes(ymin=AME-1.96*SE,ymax=AME+1.96*SE),width=0,size=.8,position = position_dodge(width=.25)) +
  theme_bw() +
  labs(x="Former vote choice",y="Effect on PTV") +
  ggtitle(" ") +
  ylim(c(-.09,.25)) +
  scale_color_manual(values=c("black","gray")) +
  theme(legend.position = "none")

ggsave("fig3r.pdf",width=4,height=4)

## Fig 3b -----

aiintm_biv_lb <- lm(maxlbptv~antiimm*expcondition,data=acdf)
aiintm_adj_lb <- lm(maxlbptv~antiimm*expcondition+gender+agegrp+education+partygroup15,data=acdf)

pgintm_biv_lb <- lm(maxlbptv~expcondition*partygroup15,data=acdf)
pgintm_adj_lb <- lm(maxlbptv~expcondition*partygroup15+gender+agegrp+education,data=acdf)

aiintmdf_lb <- bind_rows(getintmarginsdf(aiintm_biv_lb,moderator="antiimm",refadv = F,modeltype="Bivariate",comptype="Accomodative v. Control"),
                      getintmarginsdf(aiintm_adj_lb,moderator="antiimm",refadv = F,modeltype="Adjusted",comptype="Accomodative v. Control")) %>% 
  mutate(model=fct_relevel(model,"Bivariate"),
         comp=fct_rev(comp))

pgintmdf_lb <- bind_rows(getintmarginsdf(pgintm_biv_lb,moderator = "partygroup15",refadv = F,modeltype="Bivariate",comptype="Accomodative v. Control"),
                      getintmarginsdf(pgintm_adj_lb,moderator = "partygroup15",refadv = F,modeltype="Adjusted",comptype="Accomodative v. Control")) %>% 
  mutate(model=fct_relevel(model,"Bivariate"),
         pg=fct_relevel(partygroup15,c("Left bloc","Other right")),
         comp=fct_rev(comp)) %>% 
  filter(pg %in% c("Left bloc","Other right","DPP"))

#get plot yaxis limits
min(aiintmdf_lb$lower)
max(pgintmdf_lb$upper)

#plots
ggplot(data=aiintmdf_lb,aes(x=antiimm,y=AME,group=model,color=model)) +
  geom_hline(yintercept = 0,linetype="dotted") +
  geom_point(position = position_dodge(width=.05)) +
  geom_rug(data=acdf,aes(x=antiimm,y=0),sides="b",alpha=.04,position="jitter",inherit.aes = F) +
  geom_errorbar(aes(ymin=AME-1.65*SE,ymax=AME+1.65*SE),width=0,size=1.2,position = position_dodge(width=.05)) +
  geom_errorbar(aes(ymin=AME-1.96*SE,ymax=AME+1.96*SE),width=0,size=.8,position = position_dodge(width=.05)) +
  theme_bw() +
  labs(x="Anti-immigration attitude",y="Effect on PTV") +
  ggtitle("Highest left bloc propensity to vote") +
  ylim(c(-.095,.265)) +
  scale_color_manual(values=c("black","gray")) +
  theme(legend.position = "none")

ggsave("fig3bl.pdf",width=6,height=4)

ggplot(pgintmdf_lb,aes(x=pg,y=AME,group=model,color=model)) +
  geom_hline(yintercept = 0,linetype="dotted") +
  geom_point(position = position_dodge(width=.25)) +
  geom_errorbar(aes(ymin=AME-1.65*SE,ymax=AME+1.65*SE),width=0,size=1.2,position = position_dodge(width=.25)) +
  geom_errorbar(aes(ymin=AME-1.96*SE,ymax=AME+1.96*SE),width=0,size=.8,position = position_dodge(width=.25)) +
  theme_bw() +
  labs(x="Former vote choice",y="Effect on PTV") +
  ggtitle(" ") +
  ylim(c(-.095,.265)) +
  scale_color_manual(values=c("black","gray")) +
  theme(legend.position = "none")

ggsave("fig3br.pdf",width=4,height=4)

